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A systematic approach is presented for developing entropy stable (SS) formulations of any order for the 
Navier-Stokes equations. These SS formulations discretely conserve mass, momentum, energy and satisfy a 
mathematical entropy inequality. They are valid for smooth as well as discontinuous flows provided sufficient 
dissipation is added at shocks and discontinuities. Entropy stable formulations exist for all diagonal norm, 
summation-by-parts (SBP) operators, including all centered finite-difference operators, Legendre collocation 
finite-element operators, and certain finite-volume operators. Examples are presented using various entropy 
stable formulations that demonstrate the current state-of-the-art of these schemes. 

I. Introduction 

A current organizational research goal of the Revolutionary Computational AeroSciences sub-project (Aeronautical 
Sciences Project) is to develop next generation high-order numerical algorithms for use in large eddy simulations (LES) 
and hybrid Reynolds-averaged Navier-Stokes (RANS)-LES simulations of complex separated flow. These algorithms 
must be suitable for simulations of highly nonlinear turbulence models across the subsonic, transonic and supersonic 
speed regimes. 

Most high-order techniques experience a loss of robustness when the solution contains discontinuities or even 
under-resolved physical features. Although a variety of mathematically rigorous stabilization techniques have been 
developed for second-order methods (e.g. total variation diminishing (TVD) limiters, * 1 and entropy stability 2 ), ex- 
tending these techniques to high-order formulations has been problematic. A typical consequence is loss of design 
order accuracy at local extrema or insufficient stabilization. It is possible to achieve high-order design accuracy away 
from captured discontinuities, and maintain sharp “nearly monotone” captured shocks. The schemes that deliver these 
features are the essentially nonoscillatory (ENO) 3 - 4 and weighted ENO (WENO) 5 ' 6 schemes. Unfortunately, nonoscil- 
latory schemes experience instabilities in less than ideal circumstances (i.e., curvilinear mapped grids or expansion of 
flows into vacuum). Because these schemes are largely based on stencil biasing heuristics rather than mathematical 
stability proofs and the theory that does exist is not sharp, 7 - 8 there is little to guide further development efforts focused 
on alleviating the instabilities; that is, until recently. A general procedure for developing entropy conservative and 
entropy stable schemes of any order appears in reference 9, and is applicable for broad classes of spatial discretization 
operators. 

Entropy stability guarantees that the thermodynamic entropy is bounded for all time in Li, provided that density 
and temperature remain positive and boundary data is well-posed and preserves the entropy estimate. Nearly three 
decades ago, entropy conservative schemes that discretely satisfy an entropy conservation property were constructed 
by Tadmor 2 - 10 for second-order finite- volume methods. These schemes were extended to high-order periodic domains 
by LeFloch and Rhode. 11 Finding a computationally efficient discrete entropy flux was a major obstacle, that was 
alleviated recently for the Navier-Stokes equations through the work of Ismail and Roe. 12 A methodology for con- 
structing entropy stable schemes satisfying a cell entropy inequality and capable of simulating flows with shocks in 
periodic domains was developed by Fjordholm et al. 13 Recently, Fisher and Carpenter 9 present multi-domain proofs 
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of entropy conservation and stability based on diagonal norm summation-by-parts (SBP) operators, that yield en- 
tropy stable methods on finite domains. Generalization to arbitrary Cartesian domains follows immediately using 
simultaneous -approximation-term (SAT) penalty type interface conditions 14 between adjoining domains. 

A remarkable property of the work in reference 9 is the generality of the entropy stability proofs. Although the 
primary focus therein is entropy stability of WENO finite -difference schemes, all proofs generalize to a broad class of 
SBP-SAT operators applied to the conservation form of the equations. Indeed, any spatial discretization that can be 
expressed as a non-dissipative, diagonal norm, SBP-SAT operator can be implemented in an entropy conservative and 
stable fashion. Thus, now it is possible to construct entropy conservative and stable formulations of arbitrary order 
for many popular discrete operators; centered finite-differences and Legendre spectral collocation being important 
examples. 

Entropy analysis of the Navier-Stokes equations appears in the early work of Hughes et al. 15 in the context of 
Galerkin and Petrov-Galerkin finite element methods. Entropy stability follows immediately if the equations are 
rotated into symmetric (nonconservative) form and discretized using the FEM. Finite-element entropy stability proofs 
for nonlinear equations written in conservation form, appear in many texts (e.g., see Hesthaven and Warburton 16 for a 
discussion of Burgers’ equation), and are valid for elements of arbitrary order. Extension to the compressible Navier- 
Stokes equations in conservation form, however, has been difficult to achieve. Indeed, a fundamental obstacle in 
finite-element method (FEM) proofs is the requirement for integral exactness, a property that is all but impossible to 
satisfy for the compressible equations. * 2 Entropy stability proofs that are based on diagonal norm, SBP-SAT schemes, 
however, do not suffer this limitation. 

The focus herein is a “high-level” overview of the mathematical concepts of entropy stability for diagonal norm, 
SBP-SAT operators, including illustrative examples that support the theory. No attempt is made to develop new 
concepts or present new SBP-SAT operators or entropy stability theory. The pivotal works (theorems) that facilitated 
recent contributions are presented, supplemented with the citations necessary to aide further investigation by interested 
readers. Only discretizations consistent with the Lax-Wendroff theorem are considered. In section II the Navier-Stokes 
equations are introduced followed by a comprehensive introduction to diagonal norm SBP-SAT theory and operators. 
The section concludes with a survey of important discrete operators that can be expressed as diagonal norm, SBP 
operators. Section III introduces the concepts of entropy consistency and entropy stability at both the continuous and 
semi-discrete level. Details of the implementation of entropy stable schemes are presented in section IV, including a 
description of the comparative technique used to build an entropy stabilized conventional algorithm (e.g., SSWENO). 
Results from several different discrete entropy stable operators are presented in section V followed by conclusions in 
section VI. Appendices provide additional implementation details. 

II. Theory: Diagonal Norm SBP-SAT Operators 

A. Governing Equations 

Consider the calorically perfect Navier-Stokes equations expressed in the form 

* + (/% = (/ (v) %, *en, t e 

Bq = gb, x e dQ. t e [0,°°), (H-l) 

q(x,0) = go(x), xeQ, 

where the Cartesian coordinates, x = (x i,X2, X3 )' , and time, f, are independent variables, and index sums are implied. 
The vectors q, /', and f <vl ‘ are the conserved variables, and the conserved inviscid and viscous fluxes, respectively. 
Without loss of generality, a three dimensional box 

^ = 4.4] x [4 a?] x 4.4] 

is chosen as our computational domain with d£2 representing the boundary of the domain. The boundary vector 
gb is assumed to contain well-posed Dirichlet/Neumann data. We have omitted a detailed description of the three- 
dimensional Navier-Stokes equations, that may be found elsewhere. 17 

a Recasting the equations in entropy variables or using non-conservation forms of the equations are not viable options when simulating high- 
speed flows because the Lax-Wendroff (L-W) theorem is not applicable in these forms. 
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B. Summation-by-parts Operators 


First derivative operators that satisfy the summation-by-parts (SBP) convention, discretely mimic the integration-by- 
parts property 

XR *R 

j <| )q x dx= §q\ x x R L - J § x qdx, (II.2) 

X L X L 

with (]) an arbitrary test function. This mimetic property is achieved by constructing the first derivative approximation, 
£>(]), with an operator in the form 


!D = fp- 1 Q,, P=P T , 0, Ct^O, 

Qj = ‘B-Q, H = diag(— 1,0, ... ,0, 1) , 


(II.3) 


with C an arbitrary vector. The matrix P can be thought of as a mass matrix (or integrator) much like in FEM, or 
a volume that contains local grid information in the context of finite-volume or finite-differences. The nearly skew- 
symmetric matrix Q, is an undivided differencing operator; all rows sum to zero, as do all columns save the first 
and last which sum to —1 and 1, respectively. The special structure of (^guarantees conservation as is proven in the 
following lemma. (The proof to this lemma can be found elsewhere. 18 ) 

Lemma 1. All differentiation matrices, ‘D, satisfying the SBP convention given in equation II. 3 are discretely conser- 
vative in the P-norm. 


While the matrix P need not be diagonal, only diagonal norm SBP operators are considered herein; only they can 
be manipulated into the entropy stable SBP-SAT (SS-SBP-SAT) form. 

The accuracy of the first derivative operator i), can be expressed as 

<tfc(x) = PXSf + Pp+x, (II. 4) 

where P p +\ is the truncation error of the approximation. Integration in the approximation space is conducted using an 
inner product with the appropriate integration weights contained in the norm P, 

XR 

J §q x dx~§ r P£>q, <() r = ((Kxi),(Kx2),...,<KxA0) r , Q = (#(*i)>?(* 2 )> • ■ • a{xn)Y , (H.5) 

XL 

with <() and q the projections of continuous variables onto a grid r ] , ■ ■ • . r,v (fully defined in the next subsection). 
Substituting equation II. 3 into II. 5, the mimetic SBP property is demonstrated, 

^PP- 1 Qci = 4> r (<B - Qj) q = <Mjv - 4>i9i ~ <\> T, D T V<1- (H.6) 


1. Complementary Grids 

Define on the interval —1 < x < 1, the vectors of discrete solution points 

X = ,x N -i,x N ] T ; -1 < xi,x 2 ,--- ,x N -i,x N < 1 (H-7) 

Since the approximate solution is constructed at these points, they are referred to as solution points. It is useful to 
define a set of intermediate points prescribing bounding control volumes about each solution point. These (N + 1) 
points are referred to as flux points as they are similar in nature to the control volume edges employed in the finite- 
volume method. The distribution of the flux points depends on the discretization operator. The spacing between the 
flux points is implicitly contained in the norm P\ the diagonal elements of P are equal to the spacing between flux 
points, 

X = (x 0 ,Xl, . . ,X N ) , X() = X I , Xn=Xn, 

(ll.o ) 

Xi - Xi- 1 = fp(, •)(,■) , i = 1,2, ... ,IV. 

In operator notation, this is equivalent to 

Ax = !P1 ; l = (l,l,...,l) r , (II.9) 

and A is as defined in equation II. 12. Note that in equation II. 8, the solution and flux coincide at the first and last points 
and thus 

/o = /(<?! ), In = fidN ) ■ (II. 10) 

This duality is needed to define unique operators and is important in proving entropy stability. 
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2. Telescopic Flux Form 

All SBP derivative operators T> can be manipulated into the telescopic flux form. 


fx( q) 

= T~ 

'Qf+Vi 

: P 

‘Af+Tp+i 

where the N x (N + 1 ) matrix A is defined as 
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that calculates the undivided difference of the two adjacent fluxes. The existence of a telescopic form for all SBP 
operators is reiterated in the following lemma, presented without proof. (The original proof appears elsewhere. 18 ) 

Lemma 2. All differentiation matrices that satisfy the SBP convention given in equation II. 3 are telescoping operators 
in the norm fP and can be expressed as in equation 11.11. 

This telescopic flux form admits a generalized SBP property. All SBP operators defined in equation II. 3 can be 
manipulated to transfer the action of the discrete derivative onto a test function with an equivalent order of approxi- 
mation. The telescopic flux form defined in equation II. 11 combined with the flux consistency condition results in a 
more generalized relation, 


Af = § T (fB- A)f = f(q N )ty N - f{q i Hi - <|>Af, (11.13) 

where 
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gVA = <tt T x +0{N- 1 ) 


with bx the local grid spacing. This is equivalent to the commonly used explanation of summation-by-parts in indicial 


form, 


N 


N - 1 


= f{qN)§N-f{q i)<t>i- 

i=l i=l 


(11.14) 


The action of the derivative is still moved onto the test function but at first order accuracy. Note that although this 
generalized property is used herein to construct entropy conservative fluxes, it is also instrumental for satisfying the 
Lax-Wendroff theorem 19 in weak form. 


3. Variable Coefficient Second Derivative Approximation 

The viscous approximations for regularized conservation laws, written in general as 

mx)v x (x)) x = (11.15) 

must also satisfy the SBP condition. (The discrete second derivative Thlffit) operators on the vector v.) Integration by 
parts yields 

XR *R 

J cj) {$v x ) x dr = <j)dv;t|^ - J (Ih'di'xdr. (11.16) 

XL X L 
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It is trivial to show that two applications of the first derivative operator satisfy the SBP condition. This is the preferred 
method for forming the viscous terms in spectral collocation. The double application of the first derivative is not 
advisable when dealing with finite difference or finite-volume operators, as the approximation using two first derivative 
operations requires a much wider stencil, is less accurate, and only leads to neutrally stable approximations. 20-22 A 
finite-difference narrow stencil viscous operator is defined as 

‘D 2 (V)=P~ l (■&) + $['&]£>), = M(p) T , [$] =diag(>(x)), 

fMm> o, fm>o, k <1U7) 

Either approach can be used to discretize equation II. 1 and leads to an expression of the form 

(•&q x (x)) x « p-‘ (-£> r rP[tf]£> + «[$]£>) q = (11.18) 

Note that this form satisfies a telescoping conservation property that is identical to that of the inviscid terms. 

4. The Semi-Discrete Operator 

Based on the previous discussion of SBP operators and their equivalent telescoping form, the semi-discrete form of 
equation II. 1 becomes 


q, = -®,[f'(q)] + 'D i [c] i j'Djq + ¥~ l g (-? +f (v) ') + V~ l g b , 
q(x,0)=goW, xGfl, 

with g b containing the enforcement of boundary conditions. Full implementation details, including the viscous Jaco- 
bian [c] ■ ■ tensors are included in previous work 18 and elsewhere. 23-27 

Remark. It is not necessary to implement an SBP scheme in flux form, but this is the natural form to add dissi- 
pation while retaining consistency with the Lax Wendroff theorem. 18 Furthermore, the semi -discrete entropy analysis 
presented in Section B relies on the existence of the flux form. 

C. SAT Penalty Boundary and Interface Conditions 

The method of imposing boundary data is critical in all numerical methods. The manner in which these conditions are 
imposed greatly affects the stability and accuracy of solutions. Accurate, stable, and conservative interface coupling 
techniques are also essential in a multi-domain (element) setting. 

A straightforward method that permits formal analysis and maintains design-order accuracy is the Simultaneous 
Approximation Term (SAT) penalty method. 2 ~ 1-27 The approximate spatial integration of the semi-discretization in 
equation 11.19, 

jl T Vq = f 0 -f N + f { ; ] - /d v) + 1 T (g b + g ,) , (11.20) 

illustrates the purpose of the penalty (g/ are the SAT interface penalty terms), that may be thought of as a technique 
for replacing some of the computed data in the approximation with known data from the boundary condition. This 
technique is mathematically well-posed for all SBP operators. 

D. SBP Operators of Different Flavors 

The foundational concepts of summation-by-parts spatial operators were introduced in the landmark paper by Kreiss 
and Scherer. 28 Since then, all three commonly used discretization approaches: finite-difference, finite-element, finite- 
volume, have developed their version of a diagonal norm SBP operator. 

1. Finite-Difference: Centered Operators of Order 2 P . 

The first to take full advantage of the SBP formalism was the finite-difference community. High-order diagonal norm, 
central difference operators with formal SBP boundary closures were first derived by Strand. 29 For example, the 
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classical operator that is fourth-order in the interior of the domain and is closed at the boundaries with second-order 
operators (i.e., ©2-4-2) can ex P resse d in matrix form as b 
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_3_ ± 59 24 
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(H.21) 


The skew-symmetric matrix Q, follows immediately from the manipulation T ‘1). 

Countless perturbations to this operator have been developed. Examples include optimized boundary closure 
coefficients, 30 and boundary closures for the dispersion-relation-preserving (DRP) interior operator. 11 

A noteworthy extension that significantly contributes to the development of entropy stable SBP operators appears in 
the works of Yamaleev et al. 8, 18,32-34 A general strategy is presented to construct Energy Stable Weighted Essentially 
Non-Oscillatory (ESWENO) finite-difference schemes, including boundary closures that maintain, wherever possible, 
the WENO stencil biasing properties and satisfy the SBP operator convention. Stability is guaranteed in the fP norm. 
The ESWENO schemes are constructed using the following four steps 

1 . Develop a finite-domain target scheme that is stable, conservative, and accurate for smooth flows by using the 
SBP framework. 

2. Recast the target scheme in the “dual grid” conservative framework of the conventional WENO approach: the 
solution is stored and advanced at the grid points while the interface fluxes are constructed at “half points”. A 
special set of flux points and interpolants are required near the boundaries to accomplish this task. 

3. Develop a finite-domain WENO biasing strategy that allows all stencil weights to deviate from their target 
values. Precise control of the the biasing mechanics ensures design-order accuracy for smooth solutions and 
essentially non-oscillatory properties at discontinuities. 

4. Test and stabilize the scheme by using a design-order, nonlinear damping term that ensures linear L? -stability of 
the WENO operator. 

A key contribution from Yamaleev’s work 8 is the recognition that a novel set of nonuniform flux interpolation points 
is necessary near the boundaries to simultaneously achieve: 1) accuracy, 2) the SBP convention, and 3) WENO stencil 
biasing mechanics. As shown in section B it is the existence of the flux-point representation of the SBP operator that 
enables the proof of entropy stability. Details about the WENO operator in general and specifically entropy stable 
WENO operators is given in appendix B. 


2. Finite-Element: Spectral Collocation 

Summation-by-parts operators also exist for nonuniform point distributions. Reference 35 proves that Legendre spec- 
tral collocation schemes can be implemented on any distribution of points in the element, provided that appropriate 
boundary conditions (or interface conditions) are implemented to guarantee conservation and stability. A brief sum- 
mary of the mechanics of diagonal norm spectral collocation operators 35 is included herein. 
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Figure 1. The one-dimensional discretization for p = 4 Legendre collocation is illustrated. Solution points are denoted by • and flux points 
are denoted by x . 


Consider numerical methods collocated at the Legendre Gauss-Labotto (LGL) points, which include the end points 
of the interval. The complete discretization operator for the p = 4 element is illustrated in Figure 1. Define the 
Lagrange polynomials relative to the discrete points, x, as 


Li(x) 


L(x) = (x — X 2 )(x — X 3 )...(x 

(l--v)L(.v) . j / \ (1+-v)L(.y) . j / x 

2L(-i) 2/i+n 


2L(+1) 


-X„-2)(x-X N -l) 

(x-Xj)L'(xj) 


2<j<N-l 


( 11 . 22 ) 


Assume that a smooth and (infinitely) differentiable function fix) is defined on the interval —1 < x < 1. Reading 
the function / and derivative f at the discrete points, x, yields the vectors 

f ( x ) = [/(*i),/(x 2 ),--- J{xn-i)J(x n )] T ; f(x) = \f(xi),f(x 2 ),--- ,f'(x N -i)J'(x N )] T (II.23) 

The interpolation polynomial /n(x) that collocates f(x) at the points, x, is given by the contraction 

fix) m f n {x) = [L(x)] r f(x) (11.24) 

Theorem 3. The derivative operator that exactly differentiates an arbitrary p ,h -order polynomial at the collocation 
points, x, is 

V = [Lj(xj)\ (11.25) 

The proof appears in any common text on spectral methods. An equivalent representation of the differentiation 
operator can also be used, which satisfies all the requirements for being an SBP operator. 

Theorem 4. The derivative operator that exactly differentiates an arbitrary p ,h order polynomial (p = N — 1 ) at the 
collocation points x can be expressed as 

D = <£- x iff (11.26) 

Only an outline of the proof is presented. First note that in addition to the equation 11.25, the exact derivative ‘-^1 
of the function /(x) may be approximated by 

fix) » = [L(x)] r f'(x) (11.27) 

A Galerkin statement demands that the integral error between the two expressions be orthogonal to the basis set, which 
in this case are the Lagrange polynomials, L(x). This statement may be expressed as 

1 

J L(x) ([L(x)] r f'(x) - [L(x)] r f(x)) dx = 0 (11.28) 


b The nomenclature (2-4-2) signifies that boundaries/interior stencils are second- and fourth-order accurate, respectively. The resulting operator 
is globally third-order accurate in L 2 . 
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or equivalently 


V? (x) = Q,f(x) 

with 

:P = f\h{x)[L,{x)] T dx 

Q = f-iL(x)[L'(x)] T dx 

Equation (11.26) follows immediately when T is symmetric positive definite (SPD), and therefore invertible. (The 
proofs that fP SPD and Q | Qj = B appear elsewhere. 35 ) 

A diagonal norm Legendre collocation operator can be constructed by approximating the integrals in equations 
11.30 by the LGL quadrature formula. Let T| = (t|i ,1)?, • ■ ■ r| jV -i ■ q y ) be the nodes of the LGL quadrature formula (i.e., 
the zeroes of the polynomial p' nl (x){\ — x 2 )), 36 and let ft)/, 1 < l < N be the quadrature weights. 

The mass and stiffness matrices, fP and Q,., are defined by the expressions 

V = I/L(r|/;x)[L(r| / ;x)] 7 'co,fifx 

Q. = IL/L(ri/;x)[L (ri/;x)] to \dx 


( 11 . 29 ) 

( 11 . 30 ) 


The matrix fP is SPD for any x. 35 

Theorem 5. The matrix T is diagonal for the case when the collocation points are the LGL quadrature points, i.e., 
x = T|. Furthermore the diagonal coefficients of B are the integration weights, CO/, 1 < l < N used in the quadrature. 

The proofs for all theorems presented in this section appear elsewhere. 35 The discrete operators are provided for 
polynomial orders one through four in appendix A. 

3. Finite-Volume: Cell Vertex Methods 

Linite volume operators are not the focus of this work. Note, however, that SPB operators for the commonly used edge- 
based finite volume approximation of the Laplacian is developed in the work of Svard and Nordstrom 37 ' 38 in which 
accurate and stable boundary conditions are constructed for general unstructured grids. The boundary conditions are 
imposed weakly in a stable and accurate manner, using a penalty formulation. The approach is valid for general 
unstructured grids. 

E. Split-form Operators: The Remarkable Q, Matrix 

The breadth of section B (defining and describing the SBP operators), is testament to the remarkable properties em- 
bedded within the matrix operators fP and Cf The (f matrix has another property that is instrumental in developing 
the general proofs of entropy stability. A brief summary of “split-form” operators that appears elsewhere 18, 39 is now 
presented. 

Consider the Navier-Stokes equations. The nonlinear product in the convective terms can be manipulated using 
the chain rule into many forms: 1) conservative (a.k.a. divergence), 2) primitive, 3) skew-symmetric, as well as others, 
with each form exhibiting its own semi-discrete characteristics (e.g. conservation, accuracy, and nonlinear stability). 
Some splittings of nonlinear conservation laws (e.g. the splitting proposed by Honein and Moin 40 ) deliver enhanced 
discrete robustness. 

Lax and Wendroff 41 rigorously established that simulations of conservation law equations must be performed using 
the conservative form of the equation with a conservative discrete operator to approximate a weak solution. It would 
appear at first glance that split-form discrete operators are not discretely conservative because they are not derived 
strictly from the divergence form of the continuous equations. Remarkably, this is not the case for discrete operators 
constructed from diagonal norm SBP operators. 

References 18 and 39 prove that any linear splitting of the divergence and product-rule forms of nonlinear conser- 
vation laws, if discretized using any diagonal norm, summation-by-parts (SBP) operator, can be cast as a telescoping 
operator consistent with the divergence form of the conservation law. As such, any of these combinations is suitable 
for simulation of the conservation law, even in the presence of shocks. 

The split-form conservation property is demonstrated in a simple example. Consider the nonlinear equation 

du I df(u ) _ Q > n f)<r<1 

dt + dx ~ u ’ - u > U - X -T (11.32) 

m(x, 0) = uo(x ) , u(0,t) = bco(t) , u(\,t) = bc\(t ) 
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where u = u(xf ) is the continuous solution vector, / = /(h) is a nonlinear flux vector of the form f(u) = v(u)w(u). 
Because the flux / in equation 11.32 is the product of two functions, the divergence term , may be approximated 
using a split operator. 

u, + af(u) x + (1 — a )(v(u)w(u) x +w(u)v(u) x ) = 0. (11.33) 

The resulting semi -discrete equation is 

Ur + a‘D{ c Wy) + (1 - a) (VDw + ‘WDv) = 0, 

v= (v(ui),v(u 2 ),...,v{u n )) t , “V = diag(v), (11.34) 

w = (w(ui),w(u 2 ),.--,w(u n )) t , TF = diag(w). 

Although the split operator is a combination of a conservation part, ©(Tf'v), and a product rule part, ( /'i)w + W'Dw), 
no explicit smoothness constraints are placed on /. 

With these definitions (and because the diagonal norm fP commutes with ‘V and ‘W), the discrete approximation 
of equation 11.34 takes the form 

u, + atP _1 Qf+ (1 — a)tP _1 ('F'Qw-l- ‘WQy) = 0 . (11.35) 


Reference 39 proves that all split-form SBP operators have the following properties. 

Theorem 6. The discrete split-form operator (see equation 11.34) can be manipulated into the telescoping form 

aP ~ 1 QfWy + ( 1 — a) ( 'F'tP - 1 Qw + TVtp- 1 Qv) = 1 Af 

for any diagonal norm SBP operator that can be expressed in the form of equation II. 3 and for any value of the 
parameter, a. 

The first term is already in conservation form and therefore satisfies the telescoping form. Manipulating the chain 
rule term produces the following expression. 


tp- 1 (‘J/p- 1 Qw+‘M/P~ 1 Qv) = fp-'A 


JV-1 r(i) 

+ £ Y, qtj+t f i,i+l V’ 1 

1 = 1 1=1 


= tP -1 Af, 


(11.36) 


with the coefficient qi,j+t\ corresponds to the (i, i + £) row and column in (/, respectively. Note that the bracketed 
term is the manipulated flux f c . Further manipulation of equation 11.36 produces a more insightful expression for the 
telescoping flux 

r k 

fj = LL q (j+l-k,j+l ) ( w i+l v j+l—k + W j+l-kV j+1 ) , 

k=u=i (11.37) 

!<; + / j + l-k<N, l<j<N-l, 

fo = WlVl, fN = W N V N . 


Reference 39 extends the analysis further by showing that the split flux form given in equation 11.37 is consistent with 
the original conservative flux and has compact support. Thus, the combined term a QWv+ (1 — a) ( ‘V Qw + ‘W(fy ) 
can be expressed as a telescoping conservative flux, has compact support (if the Q, is compact), and is consistent 
with the original conservative flux in equation 11.32. All of the sufficient conditions of the L-W Theorem are met, so 
converged solutions using the above split operators are then weak solutions to the conservation law. 

Remark. The conservative flux fj in equation 11.37 is composed of the weighted sum of local dyadic fluxes 
(wj+iVj+i-k + w j+i-kV j+i), a form closely related the entropy conservative flux introduced later in the section C. 


III. Entropy Consistent and Entropy Stable SBP Operators 

A. Continuous Analysis 

Consider a nonlinear system of equations (e.g. the Navier-Stokes given in equation II. 1), and assume that the solution 
is smooth for all time. The objective is to bound the solution as sharply as possible. A quadratic or otherwise convex 
extension of the original equations is sought, that when integrated over the domain only depends on boundary data 
and dissipative terms. Fortunately, the Navier-Stokes equations have a convex extension, referred to as the entropy 
function that provides a mechanism for proving stability of the nonlinear system. 
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Definition 1. A scalar function S = S(q) is an entropy function of equation 11.1 if it satisfies the following conditions: 

• The function S(q) is convex and when differentiated, simultaneously contracts all spatial fluxes as follows 

V*, = S q f q q Xi = F‘q Xi =F‘_ ; i = 1 , • • • , d (III. 1 ) 

for each spatial coordinate, d. The components of the contracting vector, Sq, are the entropy variables denoted 
as w T = Sq. F‘(q) are the entropy fluxes in the i-direction. 

• The entropy variables, w, symmetrize equation 11.1 if w assumes the role of a new dependent variable (i.e., 
q = q(w)). Expressing equation 11.1 in terms of w is 

% + ( /% - (f M % = q w w t + (. fl)w Xi - (< dijW Xj ) xj =0 ; i = 1, • • • ,d (III.2) 

with the symmetry conditions: q w = [ q iV \ T , f l w = f l J ,c, ; - = c' 1 ,]. 

Because the entropy is convex, the Hessian Sqq = w q is symmetric positive definite, 

o, VC^O, (IIL3) 

and yields a one-to-one mapping from conservation variables, q, to entropy variables, w T = S q . Likewise, w q is SPD 
because q w = Wq~ l and SPD matrices are invertible. The entropy and corresponding entropy flux are often denoted 
an entropy-entropy flux pair, ( S,F ). Likewise, the potential and the corresponding potential flux (defined next) are 
denoted a potential-potential flux pair, (cp.X) t). 10 

The symmetry of the matrices q w and f‘ w , indicates that the conservation variables, q, and fluxes, /', are Jacobians 
of scalar functions with respect to the entropy variables, 

/ = < pw, [f] T = Y W , (in .4) 

where the nonlinear function, cp, is called the potential and t f are called the potential fluxes. 10 Just as the entropy 
function is convex with respect to the conservative variables ( Sqq is positive definite), the potential function is convex 
with respect to the entropy variables. 

The two elements of Definition 1 are closely related, as is shown by Godunov 42 and Mock. 43 Godunov proves 
that: 


Theorem 7. If equation 11.1 can be symmetrized by introducing new variables w, and q is a convex function of cp, then 
an entropy function S(q) is given by 

cp = w T q — S, (III. 5) 

and the entropy fluxes F'(q) satisfy 

\ f=w T f i -F i . (III. 6) 


Mock proves the converse to be true: 

Theorem 8. If S(q) is an entropy function of equation 11.1; then w T = Sq symmetrizes the equation. 

See reference 44 for a detailed summary of both proofs. Entropy analysis is now applied to the Navier-Stokes to 
determine the limits of nonlinear stability. 

Contracting equation II. 1 with the entropy variables results in the differential form of the entropy equation. 


S q q,+S q f{q) Xi = S t + F Xi = S q $ = (w T f {v) ) ^ 


-Kf {v)i = 






(III.7) 


Integrating equation III. 7 over the domain yields a global conservation statement for the entropy in the domain 


d 

dt . 


I Sdxj= w r /( v I — F dQ - f w Xj CijW Xi dxi. 


(HI- 8) 


It is shown elsewhere 9 that c (/ - in the last term in the integral are positive semi-definite. Note that the entropy can 
only increase in the domain based on data that convects and diffuses through the boundaries. The sign of the entropy 
change from viscous dissipation is always negative. 

Thus, the entropy equation derived in equation III. 8 is the convex extension of the original Navier-Stokes equations, 
and the entropy function serves as an estimator of system stability. 
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B. Semi-Discrete Entropy Analysis 

The semi-discrete entropy estimate is achieved by mimicking term by term the continuous estimate given in equation 
III. 8. The nonlinear analysis begins by contracting the entropy variables, w T , with the semi-discrete equation 11.19. 
(For clarity of presentation, but without loss of generality, the derivation is simplified to one spatial dimension. Tensor 
product algebra allows the results to extended directly to three-dimensions.) The resulting global equation that governs 
the semi-discrete decay of entropy is given by 

w r tPq,+w r Af = w r Af (v) +w 7 g i , (III.9) 

with 

w= (w(qi) T ,w(q 2 ) T ,...,w(q N ) T ) , 

the vector of entropy variables. Each semi-discrete term is now analyzed to demonstrate that it mimics the correspond- 
ing term in continuous entropy estimate, provided that a diagonal norm SBP operator is used. The form of the penalty 
terms is presented in a later section C. 

1. Time Derivative 

The time derivative is by definition in mimetic form for diagonal norm SBP operators. Arbitrary diagonal matrices 
commute, so the pointwise definition of entropy 


w] ( qi)t = C St), , Mi. 


yields the expression 


w r fPq, = l T Pw T q, = l T tPS u q, = l T PS t . 


2. Entropy Consistent Inviscid Fluxes 

The inviscid portion of equation III. 9 is entropy conservative if it satisfies 

w T Af=F(q N )-F{qi) = l r AF. (III. 10) 

Recall that w and f, F are defined at the solution points and flux points, respectively. One plausible solution to equation 
III. 10 is a pointwise relation between solution and flux-point data, which telescopes across the domain and produces 
the entropy fluxes at the boundaries. Tadmor 10 developed such a solution based on second-order centered operators. 
Herein, this solution is generalized to diagonal norm SBP operators of any order. 

Theorem 9. The local conditions, 

(w, +1 -Wif f) =tj/,+ i -tjh, ; tj>i=\|/i, W = W (III. 11) 

when summed, telescope across the domain and satisfy the entropy conservative condition given in equation III. 10. 
The potentials \j /, + 1 and % need not be the pointwise V |/ I+ 1 and \|/ ( -, respectively. A flux that satisfies this condition 
given in equation III. 11 is denoted f . 

Proof. Substituting the definition for generalized summation-by-parts in section 2, A = $ — A, into the global entropy 
conservation condition in equation III. 10 yields 

w r ®f- w r Af- l r $F + l r AF = w r ®f- l r $F — w r Af = 0. (III. 12) 

The boundary terms in equation III. 12 can be reorganized as 

w r ®f- l r ®F = (w T N f N - F n ) - K/i - Fi) = w - Vt = (HI. 13) 

where \|/i and \j/y represent the potential flux defined in equation III. 6. Defining [f] as a diagonal (N + 1) x (N + 1) 
matrix containing the elements of f, and substituting equations III. 12 and III. 13 into III. 10 yields 

— w r A[f])I = 0. 
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(III. 14) 


Substituting the equality tj/ 7 ®! = \|> 7 A1 into the left hand side of the equation yields 


(tj/^A — w r A[f]) 1 = 0. 


This is satisfied by the vector sufficient condition. 


\j/ r A = w r A[f], vj/i=\|/i, W = t) f N . 


(III. 15) 


A pointwise examination of the vector condition yields the desired result. 


□ 


3. Entropy Stable Viscous Terms 

Using the formalism introduced in section 3, viscous terms are now defined such that the continuous entropy properties 
are mimicked by the semi-discrete equation. This requires that the discrete viscous fluxes are written in terms of the 
discrete gradients of the entropy variables. 


The accuracy requirements are automatically satisfied. The coefficient matrices [c] and [c] are positive semi-definite 
because they are constructed using block-diagonal combinations of positive semi -definite matrices. 

The contribution of the viscous terms to the semi-discrete entropy decay rate is 


with the last two terms negative semi-definite. e Note that definite properties are not obvious if the discrete viscous 
fluxes are constructed based on gradients of primitive or conservative variables. While at the continuous level, w x = 
w q q x , at the discrete level in general we must assume THv ^ w u ‘Du . As in equation III. 8, only the boundary term will 
result in growth of the entropy, and thus this approximation of the viscous terms is entropy stable. 

Remark. Equation III. 16 is specific to a finite-difference SBP operator. The truncation term for a spectral colloca- 
tion operator is of the form t T p+ 1 , and the numerical dissipation is %.{£) = 0. 

The three semi-discrete terms in the entropy estimate, mimic the continuous estimate. 

C. Entropy Consistency of Inviscid Burgers Equation 

An entropy analysis of one -dimensional Burgers’ equation is presented before that of the Navier-Stokes equations. 
Conventional energy estimates exist for Burgers’ equation for all diagonal norm SBP operators. 16 Furthermore, an 
entropy estimate exists for the second-order diagonal norm operator. 10 The goal herein is to express Tadmor’s second- 
order entropy proof of Burgers’ equation in terms of diagonal norm SBP operators. 

Consider the inviscid Burgers’ equation given by 


with suitable initial and boundary data, discretized with a diagonal norm SBP operator. Alpha split the equation as 
described in equation 11.33 to obtain the form 


Now assign the flux splitting w = m/2, v = u to the chain rule terms, and assign the splitting parameter a = 2/3, which 


left multiplying with the norm fP and contracting the result with the discrete vector u to yield the energy equation 


{cw x ) x = T 1 A{ ( ' ) = t D 2 (c)’w + ‘Tp 2 pp, 


2T(c)w = fP 1 (— M{c) + < B[c\‘D) w, 


(III. 16) 


5W (c) = © r fP[c] © + $Jc), %{c) = £ 7\& r [c\kP<k- 



(III. 17) 


u t + (1/2 u 2 ) x = 0 


u, = -afP~ 1 QTPv-(l-a)(T'fP- 1 Qw+lVfP~ 1 Qv) = -fP _1 Af . 


(III. 18) 


is known to be the value required for a canonical splitting of Burgers’ equation. 16 Form the semi-discrete energy by 



c The last term is design-order small for finite-difference discretizations that use a narrow stencil viscous discretization 22 and vanishes for viscous 
terms constructed from two first-order operators (e.g. spectral collocation). 
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valid for any diagonal norm SBP operator. 

An alternative approach uses the mechanics of entropy analysis . 10 An entropy-entropy flux pair, a potential- 
potential flux pair, and the entropy variable, w, for Burgers’ equation are 

( IT \ ( IT \ 

(s,f) = ( Y’yJ ; = ( 2 ^’ 6^/ ’ u = w ' (in.20) 

Note that the entropy is guaranteed convex ( S uu = 1) for all u, and that the entropy is (chosen) equivalent to the energy 
used in the SBP analysis . 0 

The semi-discrete entropy estimate begins by left multiplying by the norm fP and contracting the result with the 
discrete entropy variable, u, to yield the entropy equation 

u T Pu, = -u r Af = -l T ZlAf = -l r AF = -(F n -Fi) (III.21) 


where f is given by the expression 


7 


2 1 i>(») ( “ ?+ “ i “ t+ “* ) ,i<.<«-i, 

jt= ( -+n=i 0 

7 _ 1 2 f _ 1 2 

Jo — 2 M i ’ Jn — 2 U n 


( 111 . 22 ) 


Equation III.22 follows immediately by expanding both the conservation and chain-rule operators in equation III. 18 
using equation 11.37 (and using a change of variables on the summation indices). 

Consider the tridiagonal second-order central operator D^ji_ \ for which an entropy flux is known to exist . 2 This 
implies that the vector relation, 1 T j ZlAf — AF] = 0, is satisfied. Define the dyadic flux in equation III. 22 to be 

- u 2 k +u k ui + uj 

Js(k,l) = 7 1 


and solve the vector equation, [ZlAf — AF] = 0, for the variables, F. The resulting fluxes, f and F, for the second-order 
D i- 2 -i operator are shown in Table C. Note that the entropy fluxes F in Table C can be expressed in the form 


Table 1. Pointwise fluxes, f and F, for the second-order SBP operator. 


f 

F 



2 


f s ( fl) 

"1 

3 



7(1.1) 


1 

6 

( M i + 

U\U 2 

+ u l) 

7(2,1) 

\u]_U 2 (u\ 

+ U 2 ) 

1 «i+H2 
2 

7(2,1) 

12 

1 

6 

( u 2 + 

U 2 U 3 

+ uf) 

7(3,2) 

\U 2 U?,(U 2 

+ M 3 ) 

1 «2+«3 
2 

7(3,2) 

M2+M3 

12 

1 

6 

(“3 + 

W3W4 

+ u l) 

7(4,3) 

g« 3 M 4(«3 

W 4 ) 

1 M3+M4 
+ 2 

7(4,3) 

«3+«| 

12 



“4 

2 


7(4,4) 

“4 

3 



7(4,4) 



N i j 

F ‘ = E ’L < lmi( u i + u k)fs(i,k) - -z( u l + u k)\ ;1<«<N-1, 

*=»'+! 1=1 (III. 23) 

1 3 - 1 3 

/(> — fN ~ 3 “iV 

The pointwise relation which facilitates a telescoping cancellation necessary to satisfy the vector relation [ZlAf — 
AF] = 0 is the local consistency condition 

(»i + 1 - Ui)fsi , i+ 1 = - u]) ; 1 < i<N— 1 . (III. 24) 

which is demonstrated next by a simple pointwise decomposition. 
d The entropy is not unique. 
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Consider the term w, (f,- — f, i ) in the vector relation [UAf — AF] = 0. Adding and subtracting equivalent terms 
yields the expression 

U/ (f/ 1'/ i ) — [2 \^i + 1 A Ui\ f; 2 [ M / "P Mi— 1 ] f/— 1 ] (III 25) 

- [i [«,-+! - Ml] f / + 5 [m; - Mi— l] l] • 

The first two terms on the RHS of equation III. 25 are already in telescoping form. The last two terms telescope if they 
satisfy a relationship of the form 


(«i+t ^ Ui)fsi , i+ 1 = V* “ Vi— 1 ;2 <i<N 

which for Burgers’ equation is satisfied pointwise by equation III. 24 (modulo slight changes at the end points of the 
domain). A decomposition of the fourth-order operator given in equation 11.21 yields a similar result. 

Equation III. 23 relies on an alpha-splitting of Burgers’ equation into canonical form suitable for entropy analysis. 
The Navier-Stokes equations, however, do not support a canonical decomposition based on the alpha-split flux tech- 
nique. Nevertheless, it is shown next that pointwise shuffling conditions similar in form to those used in equations 
III.24 and III. 25 are sufficient to achieve a generalized telescoping entropy flux for the Navier-Stokes equations. 

D. Entropy Consistency of the Euler Equations 

A general strategy for constructing high-order entropy conservative fluxes is presented elsewhere 9 and utilizes linear 
combinations of g/j-weighted, two-point entropy conservative fluxes. This approach follows immediately from the 
structural properties of diagonal norm SBP operators and the generalized SBP given in section 2. Because the ap- 
proach only relies on the existence of a two-point entropy conservative flux formula (e.g. equation III.24) and on the 
coefficients of the Q, it is valid for any SBP operator which satisfies the constraints given in equation II. 3. 

The proofs of this alternative approach for building entropy conservative operators of any order are quite involved. 
For brevity only two of the theorems are included herein. Interested readers should consult the reference 9 for details. 

The first theorem establishes the accuracy of the new fluxes. Specifically, that a high-order flux constructed from a 
linear combination of two-point entropy conservative fluxes, retains the design order of the original discrete operator 
for any diagonal norm SBP matrix Q. 

Theorem 10. A two-point entropy conservative flux can be extended to high order with formal boundary closures by 
using the form 

N i 

fr = E 'L 2 m ) fs(.qe,qk), 1 <«■<#- 1, an.26) 

k=i+ie=i 

when the two-point non-dissipative function from Tadmor 10 is used 

1 

fs(qk,qe) = f g(w(q k )+£,(w{qe)-w(q k ))) d£, g(w{u)) = f (u) . 

0 

The coefficient, qiyy), corresponds to the (k,£) row and column in Cf respectively. 

Proof To show the accuracy of approximation, the flux difference is expressed as 

fi S) - fi- 1 = E E 2( l((k)fs (ue,u k ) - E E 2< l(t,k)fs Of, U k ) , 2 < i < N — 1 . 

k =i+u=\ k =a=i 

that may be manipulated into the form (see reference 9) 

fi S) -t\ = t 2 q(i,j)fs(“ i ’Uj), 1 < i < N. (III. 28) 

;'= i 

This form facilitates an analysis by Taylor series at every solution point using the expression for the two-point fluxes 
given in equation III.27. The remainder of the proof is presented elsewhere. 9 □ 

The second theorem establishes that the linear combination does indeed preserve the property of entropy stability 
for any arbitrary diagonal norm SBP matrix Q. 
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Theorem 11. A two-point high-order entropy conservative flux satisfying equation III. 11 with formal boundary clo- 
sures can be constructed using equation 111.26, 

ff ] = £ i2q m f S (qe,q k ), l<i<N-l, 
k=i+u=l 

where fs(qe,q k ) is any two-point non-dissipative function that satisfies the entropy conservation condition 

(we - Wk) T fs (qe,qk) = Vr - Va- (HI-29) 

The high-order entropy conservative flux satisfies an additional local entropy conservation property, 

w r fP _1 Af (S) =‘p- 1 AF = F x (q) + ‘T d , (III.30) 

or equivalently, 

wf (/f -i?3) = (Pi-Pi- 0 » 1 <i<N, (III. 31) 

where 

N i 

F i= Y Y^ejc) [(we + w k ) T fs(qe,qk)~ (Va + Va)] . 1 < * < 2V — 1 . (III.32) 

k=i+U = l 

Proof For brevity, the proof is not included herein, but is reported elsewhere. 9 □ 

Remark. The existence of a local second-order entropy flux satisfying the two point shuffle relation given in 
equation III. 29 is a very strong constraint, and has until recently been a computation bottleneck. 12 

E. Entropy Stability of the Euler Equations 

Tadmor 10 identifies three “tools of the trade” in the analysis of entropy stability: comparison arguments, a homotropy 
approach and kinetic formulations. Herein a comparison approach is used to establish entropy stability. In a com- 
parison approach, the entropy dissipation generated by the primary scheme is compared with the baseline entropy of 
a scheme known to be at least entropy conservative. If the dissipation is less than the entropy conservative datum, 
then more dissipation is necessary. A popular example of schemes developed using a comparison approach, are the 
E-schemes of Osher 45 that use a Godunov scheme as the entropy datum. Conditions that guarantee entropy stability 
are now established. 

A condition analogous to equation III. 10 that guarantees entropy stability is 

yv T 'Pq t +F(q N )-F(q l )<w T g b , (III.33) 

which is satisfied if the “baseline” entropy stable inviscid fluxes satisfy the comparison condition 

w r Af > l r AF. (III. 34) 

Using the result for the entropy conservative flux in equation III. 10, this condition can be rewritten as 

w r Af > w r Af (S) . 

Substituting the generalized summation-by-parts property, 2 

w r («- A)(f (S) -f) = w^f-f^) = w r A([f] - [f (S, ])l < 0, (III. 35) 

yields the sufficient local condition for entropy stability 

w r A[f] < w r A[f (S) ], (III. 36) 

or in indicial form, 

(wi+\ — wi) T (f — f^) < 0, i= 1,2, ... ,1V— 1. (III. 37) 

The development of this condition relied heavily on the formalism introduced in Section II, including the generalized 
SBP property, that made it possible to extend the entropy stability condition to finite domains and multiple blocks. The 
entropy stability condition in equation III. 37 can be used to find entropy stable conditions for any type of telescopic 
flux operator. 

Remark. The entropy stability condition is based on the global property given in equation III. 34. It does not inform 
how to add sufficient dissipation such that a non-oscillatory solution is obtained. 

e The bracket nomenclature around a vector denotes a diagonal matrix with vector elements injected on the diagonal. 
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IV. Implementation Details 


A. Three-Dimensional Curvilinear Navier-Stokes 


All the analysis and proofs based on one-dimensional operators developed in sections II and III extend directly to 
three dimensions using tensor product algebra. Entropy consistency extends directly to the equations in contravariant 
form with the appropriate choice of metrics. The proofs of contravariant entropy consistency and the construction of 
metrics that satisfy the geometric conservation laws and entropy consistency, can be found in the thesis of Fisher. 17 
The semi-discretization of the three-dimensional Navier-Stokes equation written in contravariant form is 



±( v ); 3 

-\\ = L% l ^ + s b ) xeQ., 

i= 1 

Bq = g b , xGd£2, f€[0,°°), 
q(x,0)=go(x), xGft, 


t G [0,°°), 


(IV. 1) 


with 

q = [/](/ 5 ® ® % ® IQ ; D 5 , = (Is ® % ® \ 2 ® IQ 

and similar tensor product expressions for the other three-dimensional derivative operators. 


B. Entropy Consistent Euler Flux 


The inviscid terms in the discretization of the Navier-Stokes equations are calculated according to equations IV.8, 
III.26, and III.27, by using the two-point entropy conservative flux of Ismail and Roe, 12 

f£ (Ui,u i+ 1) = (p-6/, P'6/61 + 8/1/5, pVjX >2 + 8 j 2 p, P'6/63 + %j3P, P'6 jH)) T , 


D = 


a 

7T 


Vj+1 

VTi - iT 


h=R 


log( 


] 

VTi' 


p = 


Pi 

VTi 


Pi + 1 

VTm 


vtv r 


J=+ 1 : 
VTi ^ VTW 


VTjpi 

V^i+lpi+l 


s/TiPi + \Ai+lPi+l 


1 

VTi 


VTiVl 


. (vT; + vhi) iPi V^+iP'+i) 

7+1 


7-1 


log I 


Tj P, 
T i + 1 Pi+1 


)(A 


H = h + -v e v e , p = 


.VTi VTVVj 

(vTi + vhi) (VTiPi-VTi+iPi+i) 

2 (^(v^pO - 1 °g(v / ^+TP/+i)) 


(IV.2) 


This somewhat complicated explicit form is the first entropy conservative flux for the convective terms with low 
enough computational cost to be implemented in a practical simulation code. Previously, Tadmor 2 derived an entropy 
conservative flux form that required integration through phase space, but this was deemed too expensive to be practical. 


C. Multi-domain Operators 

A multi-domain framework (hexahedral elements in context of spectral collocation) greatly simplifies the grid gen- 
eration process for complex configurations by breaking the geometry into the union of piecewise smooth hexahedral 
domains. High-order SBP techniques naturally extend to multi -domain discretizations. 14 46-48 Each domain (element) 
is discretized with a stable tensor-product formulation and then connected to its adjoining neighbors by using interface 
conditions that maintain the stability, accuracy, and conservation of the interior operators. Domain interfaces only 
need be Co smooth to maintain the stability, conservation and design accuracy. 

The entropy stability proofs extend directly to a multi-domain (element) context for SAT inviscid and viscous 
interface penalties. Contracting the three-dimensional Navier-Stokes equations against the entropy variables in a multi- 
domain context, produces both inviscid and viscous boundary interface entropy flux terms in each spatial direction of 
the following form 

w r A^[f^ - fji v) ] = ^(F - F M)f . (IV. 3) 
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1. Inviscid Interface Dissipation 

Inviscid interface penalties replace ff with an entropy conservative flux V r (f Kq- ' i ) plus a dissipative stabilization 
term. The interface flux f sr is constructed to satisfy the consistency relation 

( w (-)_ vv (+))f'-(g,(-) ;? (+)) = \|/(+) — \|/(-) . 

Lax-Friedrichs dissipation is used to stabilize the interface. The resulting flux is expressed as 

fssrt q {-)^q(+h = pY^H^W) + ^fa(vv(-)-w( + )), 

I (IV 4) 

X=[i(( M H) 4 + (cW) 4 + ( M H) 4 + (cW ) 4 ) ] ? 

The dissipation term is scaled by any positive number. The resulting term that appears in the entropy estimate is of the 
form 

-MvvH-w(+)) 2 (IV.5) 


2. Viscous Interface Dissipation 

Viscous penalties are constructed following the technique presented elsewhere. 17 


D. The Primary Scheme: WENO 

A WENO scheme is used herein as the primary scheme. Note, however, that in principle any scheme may be used 
as the primary operator in a comparative approach. WENO methods are well suited for high-fidelity simulations of 
complex physics that contain discontinuities. An example of a typical application is a canonical simulation of sound 
that is generated by a shock- vortex interaction. 49 - 50 The high-order nature of WENO efficiently resolves the subtle 
details of the sound generation and propagation, while the stencil biasing mechanics ensures robust, high resolution 
properties in the vicinity of the shock. 

The WENO algorithm is quite involved and a brief description is included in appendix B. The algorithm is im- 
plemented using a dual-grid approach with uniformly spaced solution points and flux points that receive interpolated 
data. Multiple candidate fluxes are built at each flux point and then a convex combination of smooth components is 
formed. The WENO flux is calculated using the formula 


Aw) = 
J j 


,k fh 


k=\ 


. «}fj 


7 = 1 , 2 ,. 


(IV.6) 


The stencil biasing procedure ensures that stencils containing discontinuities result in high relative values of p and 
are assigned negligible weights. The WENO flux becomes an interpolation that only incorporates smooth data. Away 
from discontinuities, the flux collapses to the target flux. Near discontinuities, the flux transforms into an upwind 
operator. Note that the flux consistency condition on the first and last flux point in equation II. 10 is enforced. The 
WENO method as described has no provable stability properties. We employ a limiting procedure to ensure that it 
satisfies entropy stability. From equation III. 37, we require that 




(IV.7) 


where f( ssw ) is the entropy stable WENO flux and f' S) is the entropy conservative flux. The limiter that guarantees 
entropy stability is chosen as 


fi ssw) =;r+ffi»-fr). 5 =^+ 2 ^, 

v 7 y/b z + c- (IV. 8) 

b = (w i+l - Wi ) T (/f -fl W) ~) , c = 1(T 12 , 

where jf V ‘ is the WENO flux without correction described above. The limiting process to get the entropy stable 
property is not unique. It was chosen merely because it is smooth with respect to the solution. This same limiting 
procedure works with fluxes other than WENO. It is shown in a later section that this minor correction to the WENO 
operator has a large impact. 
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V. Results: Accuracy and Robustness 


The accuracy and robustness of the algorithms developed herein are tested using two smooth and two discontinuous 
problems. The smooth problems are the propagation of an isentropic vortex and the propagation of the viscous shock. 
Both problems demonstrate the design-order convergence of the new entropy conservative formulation and of the 
SSWENO scheme. The final two problems (the Sod and Lax shock tube problems) have discontinuous solutions, and 
test the efficacy of the entropy stable correction terms used in the SSWENO formulation. It is demonstrated that the 
corrections do not adversely affect the desirable stencil biasing properties of the baseline WENO scheme. 


A. Isentropic Vortex 

The isentropic vortex is an exact solution to the Euler equations that is an excellent test of the accuracy and function- 
ality of the inviscid components of a Navier-Stokes solver. It is fully described by 


f(x,y,z,t ) = 1— (x — xo — f/ 00 cos(a)f) 2 + (>' — yo~ ^~sin(a)f) 2 
T{x,y,z,t) = [l — £ 2 Mi^exp(/(x,y,z,f))l , P (x,y,z,t) = TF i, 


y— yo~ £/oosm(a)/ 
2n 


exp 


2 


u(x,y,z,t) = L/ooCos(a)- 

v(x,y,z,t) = t/,sin( a) ~e /-*°-^ cos(a > exp (&**> 
w(x,y,z,t) = 0. 

In this study the values U,„ = M m Coo, £ v = 5.0, ;V/ :Xi = 0.5, and y = 1 .4 are used. 


(V.l) 


1. Finite Domain Cartesian Grid Test 

The finite domain Cartesian grid test is described by 


xe (-5,5), ye (-5,5), (x 0 ,yo) = (0,0), a = 0.0, t > 0. 

Three different grid resolutions are examined, and the vortex is halfway out of the domain when the error measure is 
evaluated. This measures the effect of the boundary closure, the penalty boundary condition, and the interior scheme. 
Both the (3-6-3) and (2-4-2) entropy stable WENO finite-difference schemes are evaluated. The designed order of 
accuracy is observed in Tables 2 and 3. The error decay asymptotes towards the designed rate in each case (i.e., third- 
order and fourth-order respectively). The reduction in design order relative to the periodic test case results from the 
order p boundary closures. 

Table 2. Error convergence is shown for the (2-4-2) simulation of the finite domain isentropic vortex. 



Entropy Consistent FD 

Entropy Stable WENO FD 

Resolution 

L? error 

L 2 rate 

L°° error 

L°° rate 

Lr error 

L 2 rate 

L°° error Z 

“ rate 

32x32 

5.37e-02 

- 

3.17e-02 

- 

5.44e-02 

- 

3.49e-02 

- 

64x64 

6. 14e-03 

3.13 

3.51e-03 

3.18 

6.26e-03 

3.12 

3.51e-03 

3.31 

128 x 128 

7.43e-04 

3.05 

4.63e-04 

2.92 

7.46e-04 

3.07 

4.72e-04 

2.90 

256 x 256 

9.01e-05 

3.04 

5.95e-05 

2.96 

8.39e-05 

3.15 

5.52e-05 

3.10 


B. Viscous Shock 

This problem tests the accuracy and functionality of the viscous terms in the solver by using an exact solution to the 
Navier-Stokes equations. The nonlinear convection and viscous terms are perfectly balanced in this computation, and 
thus the shock thickness remains constant and is simply advected. The derivation of this problem is available in the 
thesis of Fisher. 17 

The convergence rate for the viscous shock is evaluated on a Cartesian grid, described by 

xe (-1,1) x (-0.5, 0.5). 
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Table 3. Error convergence is shown for the (3-6-3) simulation of the finite domain isentropic vortex. 



Entropy Consistent FD 

Entropy Stable WENO FD 

Resolution 

L? error 

L 2 rate 

L°° error 

L°° rate 

L 2 error 

L 2 rate 

L°° error L°° rate 

32x32 

7.33e-02 

- 

3.65e-02 

- 

7.32e-02 

- 

3.73e-02 

- 

64x64 

4.30e-03 

4.09 

2.10e-03 

4.12 

4.95e-03 

3.89 

3.03e-03 

3.62 

128 x 128 

1.94e-04 

4.47 

2.08e-04 

3.34 

2.53e-04 

4.29 

3.80e-04 

3.00 

256 x 256 

7.72e-06 

4.65 

1.79e-05 

3.54 

1.09e-05 

4.54 

3.09e-05 

3.62 

512x512 

2.83e-07 

4.77 

1.01e-06 

4.15 

4.19e-07 

4.70 

1.71e-06 

4.18 


The shock flow is rotated a = 20° with respect to the jc-axis. The shock profile is initially located at x s = —0.5 with 
respect to the origin, and is simulated until t = 0.25. The Reynolds number was Re = 10 and the reference Mach 
number was M = 2.5. The errors for the Cartesian grid are shown in Table 4. The errors for the (2-4-2) and (3-6-3) 
operators asymptotically converge at the designed orders of 3 and 4, respectively. 

Table 4. Error convergence is shown for (2-4-2) and (3-6-3) simulations of the viscous shock. 



Entropy Stable(2-4-2) 

Entropy Stable (3-6-3) 

Resolution 

L? error 

L 2 rate 

L°° error 

L°° rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

48 x24 

1.74e-03 

- 

5.60e-03 

- 

3.36e-04 

- 

1.57e-03 

- 

96x48 

1.10e-04 

3.97 

4.24e-04 

3.72 

7.62e-06 

5.46 

6.51e-05 

4.59 

192 x 96 

7.30e-06 

3.92 

3.38e-05 

3.65 

2.53e-07 

4.91 

3.42e-06 

4.25 


C. Shock Tube Problems 

The shock tube problems below show that no adverse effects of the entropy correction are observed for flows that 
admit shocks. 




(a) Sod Shock Tube, t = 0.2 


(b) Lax Shock Tube, r = 1.3 


Figure 2. Shock tube solutions are plotted for the entropy stable WENO methods developed in this work and compared to reference 
solutions. 
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1. Sod Shock Tube 


Sod’s shock tube problem evaluates how a numerical method behaves when a shock, expansion, and contact discon- 
tinuity are present. Of particular note is how much smearing is observed in the shock and contact, and whether any 
oscillations are present. 

The domain is 

xe ( 0 , 1 ), y e (—0.1, o.i), t> o, 

and is initialized with. 


PC x,y,z) = 


1 x < 0.5, 
1/8 x > 0.5, 


p(x,y,z) = 


y x < 0.5, 
y/ 10 x>0.5, 


u{x,y,z) = 0, v(x,y,z) = 0, w{x,y,z) = 0, 


(V.2) 


where y = 7/5. The problem is simulated with (3-6-3) and (2-4-2) entropy stable WENO operators with N = 100 
uniform cells. The solution is plotted for t = 0.2 in Figure 2(a). The solutions do not exhibit oscillations and the shock 
smearing is nearly equivalent between the two schemes, with slightly less diffusion observed in the (3-6-3) scheme. 


2. Lax Shock Tube 


Lax’s shock tube problem is used to show that no entropy problems are observed using the current methodology and 
that the correct shock location is observed. The reference solution uses N = 800 points with the (2-4-2) WENO 
operator. 

The simulated domain is 


initialized with 


xe(-5,5), y£ (-0.5, 0.5), t > 0, 


P(- x,y,z) = 


0.445 x < 0.0, 
0.5 x > 0.0, 


p(x,y,z) = 


3.528y x<0.0, 
0.571y x > 0.0, 


u(x,y,z) = 


0.698 x < 0.0, 
0.0 x > 0.0, 


v(x,y,z)=0, w{x,y,z)= 0, 


(V.3) 


where again y = 7/5. The simulation used N = 200 uniform cells and the solution is plotted in Figure 2(b) for f = 1.3. 
Again the solutions do not exhibit oscillations. 


D. Applications 

Two simulations of two-dimensional flows are presented to demonstrate the capability of the multi-block SSWENO 
algorithm. 


1. Multi-Element Airfoil 

The first problem is a simulation of a multi-element airfoil. The Mach number is M = 0.3, the angle of attack is 20°, 
and the Reynolds number is Re = 30,000 based on chord. The fourth-order SSWENO formulation is used and the 
simulation time is taken to T = 30.0. 

Although no experimental data is available for these flow conditions the complexity of the geometry tests the 
flexibility of the multi-block SSWENO formulation. The grid contains 73 point-matched blocks constructed with the 
Pointwise commercial grid generation package. 51 The small size of the blocks (e.g., O (50 x 50) points) enable locally 
analytic mesh generation in each, ensuring high-order smoothness of the mesh and thus differentiable metrics. Figures 
3(a) and 3(b) show iso-contours of the magnitude of vorticity and the entropy, respectively. Simulations with the 
SSWENO formulation are both more accurate and robust than either conventional WENO or ESWENO 8 ' 32 ' 33 based 
on the same fourth-order target operator. 
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(a) Vorticity (b) Entropy 


Figure 3. Two solutions for a multi-element airfoil are presented, demonstrating the feasibility of simulating complex flow features with 
the entropy stable WENO formulation. Red and blue vorticity contours identify regions of positive and negative vorticity. 


2. Supersonic Cylinder 

A problem adapted from Chaudhuri et al. 52 is used to simulate shock-vortex interactions in viscous supersonic flow 
around a cylinder in a duct. The geometry is simple by multi-block standards, but the bow shock can create a stability 
problem for multi-block interfaces if the treatment lacks sufficient robustness. 

The cylinder is located at (x, y)/D = (0, 0). The duct inlet is located at x/D = 3 and the duct outlet is at x/D 
= 21. The top and bottom duct walls are located along y/D = 3, respectively. An O-type multi-block configuration 
is used to discretize the domain surrounding the cylinder. Five topological blocks are needed to fully discretize the 
domain. The cylinder walls are treated as adiabatic no slip walls. The duct walls are treated as slip walls and the grid 
spacing is approximately isotropic. The inflow is a uniform freestream condition at M = 3.5. The Prandtl number 
used is Pr = 0.7 and the Reynolds number based on diameter is 10 4 . The fourth-order SSWENO algorithm is used in 
the study. A shock sensor is used to increase efficiency, by deactivating the stencil biasing mechanics in the WENO 
algorithm in regions where the solution is smooth. 

Iso-contours of the density, Mach number, entropy and the shock sensor are shown in figure 4. Note that shock- 
vortex interactions pervade the entire wake region, and the shocks also propagate and interact throughout the domain. 
Many modes of large scale unsteadiness are observed where the reflected shocks move back and forth downstream of 
the cylinder and the vortex structures propagated through the domain exhibit different pairings for different times. It 
is recognized that the true physical problem would be a three-dimensional turbulent flow. However, this simulation 
suffices for the purpose of demonstrating the high Mach number capability of the multi-block SSWENO formulation. 

VI. Conclusions 

A “high-level” overview is presented of the mathematical concepts of entropy stability for the Navier-Stokes 
equations. Recent contributions to the field are summarized that prove that all diagonal norm, SBP-SAT operators 
may be implemented in an entropy conservative (Euler) or entropy stable (Navier-Stokes) fashion. Thus, entropy 
stable operators of arbitrary order may be constructed for the Navier-Stokes equations, that guarantee an Lo bound on 
the thermodynamic entropy, provided that density and temperature remain positive and boundary data is well-posed. 

Many popular discrete operators may be represented as diagonal norm SBP operators, including all centered finite- 
difference, and Legendre spectral collocation operators. Extension to three-dimensional geometries via a multi-block 
approach with SAT interface coupling follows immediately. Some finite-volume operators also satisfy the necessary 
constraints, but are not the focus of this overview. 
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A comparative approach is used to combine an entropy conservative or entropy stable operator with a conventional 
high-order operator, thereby guaranteeing an entropy bound for the combined operator. Herein, a conventional multi- 
block WENO operator is the method of choice for simulating complex geometries with strong shocks. It is combined 
with the entropy conservative formulation to produce a multi-block entropy stable WENO operator. (In principle, any 
high-order dissipative scheme would suffice as a candidate for a comparative operator.) Important implementation 
details of the combined algorithm are included. 

Test cases demonstrate the efficacy of the multi-block, entropy stable WENO algorithm. Design order accuracy 
is achieved for smooth flows using the Euler vortex and the viscous shock test cases. Sod’s and Lax’s shock tube 
problems demonstrate the superior shock capturing capabilities of the entropy stable WENO scheme. A multi-element 
airfoil simulation demonstrates the geometric flexibility of multi-block component of the algorithm. Finally, a cylinder 
in a Mach 3.5 cross-flow demonstrates the shock capturing capabilities of the algorithm. 
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A. Spectral Collocation operators 

The discrete operators are provided for polynomial orders one through four in Table 5. Three quantities completely 
define the discrete operators. They are the diagonal norm F, the nearly skew-symmetric Cg and the positions of the 
collocation points, x. Recall the differentiation matrix is D = F 1 Q. Only the upper triangular portion of Q, is 
provided. The full Q, matrix can be reconstructed from the skew-symmetry property Q,+ Qj = B. 

B. Entropy Stable WENO Finite Differences 

Non-dissipative numerical methods cannot be used to simulate shocks. The primary reason for this is that shocks 
dissipate energy and non-dissipative numerical methods have no mechanism to mimic this. To simulate problems with 
shocks, dissipation needs to be added to the numerical method. There are a variety of mechanisms to achieve this, but 
in this work Weighted Essentially Non-Oscillatory (WENO) finite-difference methods are used. The implementation 
uses unique formal boundary closures from Fisher et al. 33 that satisfy the SBP condition. Stencil biasing mechanics 
follow two papers by Yamaleev and Carpenter. 8,32 The details of the generally applicable correction procedure are 
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Table 5. Differentiation operators for polynomials of degree one through four. 


p = 

1 

P = 

2 

P = 

3 

P = 

4 

xl = 

-1 

xl = 

-1 

xl = 

-1 

xl = 

-1 

x2 = 

+ 1 

x2 = 
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-1/V5 

x2 = 

-VV7 
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+1 
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+1/V5 
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+1 
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x5 = 
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detailed below. The full implementation details including the WENO stencil biasing algorithm throughout the domain 
are available in Fisher et al. 33 and Carpenter et al. 34 The first step to construct a WENO finite-difference operator is 
to cast the difference operator in flux form, 

Qf = Af . 

These fluxes that recover the non-dissipative first derivative approximation are called target fluxes. The target fluxes 
are broken into a sum of fluxes on smaller stencils of width, p, called candidate stencils, 

" s 

fr > ;' = i,2,...,iv-i, (B.i) 

k= 1 


where n s is the number of candidate stencils needed to describe the target flux, fj k , are the candidate fluxes, and d k are 
the target weights that recover the target flux. The candidate stencil width is held constant for all fluxes in the domain. 
For example, fourth-order operators use p = 2 and sixth-order operators use p = 3. The fourth-order case is shown in 
Figure 5. The number of candidate stencils needed to describe the fluxes, fj, can vary, as the target fluxes do not all 
have the same stencil size when approaching the boundary. The functional form of the candidate stencils depends on 
the distribution of the flux points, and thus is fully described by the norm, T. and the desired order of accuracy. 

WENO works by preventing the interpolated fluxes, fj, from using data across discontinuities. This is done by 
replacing the target weights, dp with nonlinear weights, 





k= 1, . . . ,n s . 


(B.2) 


The functional form of the nonlinear weights relies on the scaling parameter, e, and dual stencil-biasing parameters, x 
and (3. x is a measure of the smoothness over the full stencil. 


x 


i 


E 


i=l 


fd 2 P- ] u(xj) 
V dx 2 P~ l 


(■ §x) 2p 


2 


n x =n s -p. 


(B.3) 


p is a measure of the smoothness over each individual candidate stencil. 


p-\ 


Pi=E( & ) 


2 1 


k fx 

dx 1 


(*/')' 


where (pj(x) is the unique order (/? — 1) polynomial fit of the solution over the candidate stencil, /*. 
The WENO flux is calculated using the formula. 


Aw) 

J j 


" s 

k= 1 


(B.4) 


(B.5) 
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Figure 4. The density, Mach number, entropy, and shock sensor of a shock-cylinder interaction, demonstrating the shock capturing 
capabilities of the entropy stable WENO formulation. 
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Figure 5. The stencil for a WENO scheme with p - 2 and n s = 4 candidate stencils is shown. 
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